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Abstract 

Using  non-equilibrium  molecular  dynamics  for  thermal  energy  transfer,  we 
investigate  the  interfacial  thermal  conductance  between  non-covalently  inter¬ 
acting  graphene  nano-ribbons  (GNRs)  of  varying  lengths  and  widths  in  a  cross¬ 
contact  (x-shaped)  geometry.  Our  results  show  that  the  out-of-plane  conductance 
between  GNRs  can  vary  significantly  (up  to  a  factor  of  4)  depending  upon  their 
geometric  parameters.  We  observe  that  when  plotted  against  aspect  ratio,  the 
predicted  interface  thermal  conductance  values  fit  excellently  on  a  single  master- 
plot  with  a  logarithmic  scaling,  suggesting  the  importance  of  GNR  aspect  ratio 
towards  thermal  conductance.  We  propose  a  model  based  on  incorporating 
different  thermal  conductance  characteristics  of  edge  and  inner  interacting 
regions  which  predicts  the  observed  logarithmic  dependence  on  aspect  ratio.  We 
also  study  the  effect  of  graphene  edge  roughness,  temperature,  and  strain  on  out- 
of-plane  thermal  conductance  and  discuss  the  observed  results  based  on  local 
vibrational  characteristics  of  atoms  within  interacting  region,  number  of  inter¬ 
acting  phonons,  and  the  degree  to  which  they  interact  across  the  interaction 
zone. 

[S]  Online  supplementary  data  available  from  stacks.iop.org/TDM/1/025005/ 
mmedia 

Keywords:  molecular  dynamics,  graphene  nano-ribbons,  interface  thermal 
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Introduction 

Technological  advancements  in  fabrication  at  ‘nano’  length  scales  over  the  past  decade  have 
opened  up  a  gamut  of  opportunities  for  development  of  next  generation  nanoelectronics,  such  as 
nano-electromechanical  systems  [1,  2],  nano  integrated  circuits  [3],  thermoelectric  devices 
[4,  5],  thermal  rectifiers  [6],  and  others.  As  this  exciting  field  continues  to  evolve  towards 
fabrication  of  next  generation  microscale  and  nanoscale  devices,  anticipated  and  un-anticipated 
challenges  are  evolving  as  well.  Due  to  the  much  larger  surface  area  to  volume  ratio  in  such 
small  devices,  power  dissipation  at  interconnects  emerges  as  one  of  the  most  crucial  concern 
towards  determining  device  efficiency,  performance  predictability,  and  reliability.  The  energy 
transfer  through  these  contacts  occurs  via  conduction  electrons  and  thermal  vibrations,  as  well 
as  their  mutual  interaction.  This  transfer  is  often  characterized  as  electrical  contact  resistance 
(for  electrons),  interface  thermal  resistance  (for  heat)  and  electron-phonon  coupling, 
respectively.  In  an  operating  device,  these  resistances  collectively  govern  the  current  flow, 
control  temperature  distribution  within  the  device,  influence  detrimental  hot  spots,  and  thus 
dictate  its  performance. 

Owing  to  excellent  thermal  and  electrical  properties,  carbon  nanostructures  offer  a  novel 
and  cost  effective  route  towards  designing  next  generation  nano-electronic  devices.  Moreover, 
it  is  often  reported  that  these  properties  can  be  tuned  through  different  chemical  and  structural 
modifications  (such  as  functionalization,  dopants,  defects,  etc)  [7-10].  Of  particular  interest  to 
this  study  is  the  issue  of  interface  thermal  conductance  (ITC)  across  carbon  based 
nanostructures,  which  is  associated  with  the  vibrational  energy  transfer  at  the  interface  (its 
inverse  is  termed  as  interface  thermal  resistance).  In  carbon  based  electronics,  such  interfaces 
either  involve  dissimilar  material  contacts,  such  as  carbon  nanotube  or  graphene  interacting 
with  other  metals  or  semi-conductors  (also  known  as  hetero-contact)  [11,  12],  or  involve  similar 
material  contact,  such  as  CNT-CNT  cross  contact  within  CNT-mesh  in  a  CNT  based 
nanoelectronics  [13-15]. 

Within  the  framework  of  similar  material  contacts  comprised  of  carbon  nanostructures, 
several  groups,  both  experimentally  [16-21]  and  through  simulations  [22-33],  have 
investigated  thermal  interface  conductance  (or  resistance)  at  CNT-CNT  and  graphene-graphene 
interfaces  over  the  past  decade.  The  studies  offer  insights  into  thermal  transport  across  such 
interfaces,  and  how  it  can  be  tuned  by  modifying  the  interface  properties  through  alignment 
[23],  surrounding  environment  [26-28],  possible  metallization  [34],  functionalization  [32,  35], 
etc.  Nevertheless,  very  few  studies  have  explored  out-of-plane  thermal  conductance  between 
graphene  layers,  which  are  interacting  through  van  der  Waals  forces  [20,  28-31].  Most  of  the 
graphene  thermal  conductance  research  has  been  focused  either  on  in-plane  transport  of 
graphene  (such  as  effect  of  grain  boundaries  or  defects  on  in-plane  thermal  conductance  at  the 
defected  interface)  [36]  or  on  graphene  thermal  conductance  at  hetero-contacts  (such  as  a 
graphene  monolayer  supported  on  different  metallic/semiconducting  substrates)  [37-39]. 
Furthermore,  to  the  best  of  our  knowledge,  no  study  has  focused  on  exploring  interfacial 
thermal  conductance  across  single-layer  graphene  nano-ribbons  (GNRs). 

The  current  study  uses  non-equilibrium  molecular  dynamics  (NEMD)  simulations  to 
investigate  whether  intrinsic  GNR  characteristics  (length,  width,  edge  features)  and  external 
stimuli  (tensile  strain,  temperature)  can  affect  the  thermal  energy  transport  across  the  GNR 
contacts,  and  if  so,  to  what  degree.  The  understanding  of  these  out-of-plane  thermal  properties 
and  contact  energy  transport  shall  provide  a  guiding  role  towards  determining  graphene  based 
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Figure  1 .  Schematic  of  a  representative  x-shaped  GNR  pair  as  considered  in  this  study. 
L  and  W  represent  the  length  and  width  of  the  nano-ribbon,  respectively.  Color  scheme: 
fixed  carbon  atoms  (black);  hot  thermostated  carbon  atoms  (red);  cold  thermostated 
carbon  atoms  (blue);  unthermostated  carbon  atoms  (cyan). 


device  inter-connect  specifications  and  GNR  junctions  to  optimize  device  performance,  achieve 
predictable  operation,  and  reliability  in  the  use. 


Simulation  methodology 

Simulation  setup 

Figure  1  shows  a  representative  system  of  X-shaped  geometry  made  from  two  perpendicular 
GNRs.  First,  a  bi-layer  graphene  system  was  generated  using  Materials  Studio  with  planar 
dimensions  of  120xl20nm  and  an  inter-layer  spacing  of  3.4  A.  Next,  two  perpendicular 
ribbons  were  carved  out  using  an  in-house  script  with  the  desired  length  (L)  and  width  (W)  as 
shown  in  figure  1.  For  the  investigation  of  length/width  dependence  on  thermal  conductance, 
several  systems  were  generated  where  L  was  varied  from  10  to  100  nm  and  W  was  varied  from  2 
to  20  nm.  All  simulations  were  performed  using  LAMMPS  molecular  dynamics  software  [40] 
using  the  PCFF  force  field  [41].  This  force  field  has  been  successfully  employed  for  the 
investigation  of  thermal  conductance/conductivity  in  several  elemental  carbon  related  studies 
[42 — 48] .  After  an  initial  minimization,  all  systems  were  subjected  to  NVT  (canonical  ensemble) 
simulations  for  equilibration  for  200-500  ps  (depending  upon  the  system  size)  where  both  edges 
(1  nm)  at  ends  of  both  GNRs  (shown  as  black  in  figure  1)  were  fixed.  During  the  equilibration 
stage,  a  timestep  of  1  fs  was  used  for  the  simulations.  As  the  length  of  the  GNR  was  fixed,  the 
initial  C-C  bond  lengths  were  adjusted  to  make  sure  that  there  was  no  residual  tensile/ 
compressive  stress/strain  within  the  GNR  during  equilibration  by  monitoring  different  pressure 
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components1.  We  also  studied  the  effect  of  temperature  as  well  as  tensile  and  compressive  strain 
on  interfacial  thermal  conductance  for  a  representative  bi-layer  GNR  system.  Later  in  the  text,  it 
is  shown  that  such  modulation  of  external  stimuli  can  lead  to  significant  modification  of 
interfacial  heat  flow. 


Thermal  transport  simulations 


After  equilibration,  NEMD  simulations  were  used  to  calculate  thermal  conductance  across  GNR 
cross-plane  interfaces.  For  these  simulations,  the  edges  (~lnm)  on  both  ends  of  the  two 
perpendicular  GNRs  were  kept  fixed.  To  introduce  a  temperature  drop  (or  discontinuity)  at  the 
interface  interacting  through  van  der  Waals  interactions,  3  nm  sections  near  both  ends  were 
thermostated  at  either  350  K  (hot  GNR)  or  250  K  (cold  GNR)  for  all  studied  cases  as  shown  in 
figure  1.  All  simulations  were  performed  with  a  timestep  of  0.5  fs  under  the  NVE 
(microcanonical)  ensemble.  To  investigate  possible  artifacts  caused  by  the  choice  of  thermostat, 
several  thermostat  schemes  were  tested  with  heating/cooling  either  at  the  boundaries  or 
throughout  the  whole  GNR  as  well  as  using  several  thermostating  schemes  (Nose-Hoover  or 
temperature  rescaling)  with  different  adjustable  parameters  associated  with  each  scheme. 
However,  only  a  slight  difference  in  estimated  thermal  conductance  was  observed  with  no 
statistically  significant  scattering  or  trends,  as  also  reported  previously  [49].  For  all  the  results 
presented,  we  used  the  temperature  rescaling  approach  in  the  thermostated  regions  with  the 
rescaling  frequency  of  100  timesteps.  Once  steady  state  was  reached  (~500ps),  two  different 
forms  of  interface  thermal  conductance  calculated  using  the  following  set  of  equations 


A\  =  — ;  A2  = - - - 

AT  AgnrAT 


(1) 


where  Ai  (pW/K)  and  A2  (MW/m2-K)  are  total  and  area-normalized  ITC,  respectively,  Q  is  the 
heat  flow  rate,  AT  is  the  temperature  difference  between  the  two  GNRs  at  the  interface  contact 
(overlapped  region),  and  AGNR  is  the  overlapping  area  of  interaction  at  the  GNR  contact.  For 
detailed  discussion  of  the  estimation  of  this  area,  please  refer  to  supporting  information  SI.  The 
heat  flow  rate  Q  was  calculated  as  follows 


Q_ 

AAt 


1 

AAt 


(2) 


where  vpk  and  vk  are  the  velocities  of  the  thermostated  (either  hot  or  cold)  atoms  before  and  after 
the  rescaling  to  desired  temperature,  respectively  and  NB  is  the  number  of  atoms  in  the 
thermostated  region. 


3  In  order  to  achieve  zero  planar  residual  strain  for  our  systems  of  interest,  a  test  case  of  nano-ribbon  was  run  with 
different  initial  C-C  bond  lengths  (between  1.41  and  1.43)  at  300K  and  planar  pressure  was  monitored.  Based  on 
these  simulations,  the  case  with  -zero  planar  pressure  was  used  as  standard  C-C  length  for  all  studied  simulations 
presented. 
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Graphene  Nanoribbon  Length  (nm) 


(a) 


(b) 


Figure  2.  (a)  Total  thermal  conductance  (pW/K)  and  (b)  area-normalized  thermal 
conductance  (MW/m2-K)  across  different  studied  GNRs,  plotted  as  a  function  of  GNR 
length  for  different  widths.  A  power  law  scaling  of  Vi  is  shown  in  (a)  as  a  guide  to  the 
eye.  Semi-transparent  data  represent  additional  simulations  (not  discussed  in  text)  to 
further  confirmation  of  the  observed  trends. 

Results  and  discussion 

Effect  of  GNR  length  &  width 

Figure  2(a)  plots  the  total  thermal  conductance  (Ai)  as  a  function  of  GNR  length  for  different 
studied  widths.  A  few  additional  GNR  systems  were  modeled  (in  addition  to  what  is  described 
in  the  simulation  methodology)  for  a  better  estimate  of  thermal  conductivity  trends,  and  are 
shown  in  light  shaded  colors4.  A  power-law  scaling  of  Vi  is  also  plotted  as  a  guide  for  the  eye. 
As  expected,  the  figure  clearly  shows  that  for  constant  GNR  length,  the  total  thermal 

4  Additional  simulations  were  performed  to  validate  the  aspect  ratio  trend,  especially  at  the  lower  aspect  ratio,  as 
shown  in  figure  2.  For  wider  GNRs,  10  nm  and  20  nm  simulation  data  points  are  missing  due  to  system  constraints, 
e.g.  it  is  impossible  to  model  20  nm  wide  GNR  with  only  10  nm  length. 
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Figure  3.  Plot  of  the  area-normalized  thermal  conductance  as  a  function  of  aspect  ratio 
(;/)  for  different  studied  GNRs  on  a  log-log  scale.  A  power  law  scaling  of  Vi  is  shown  as 
a  guide  to  the  eye. 

conductance  increases  monotonically  with  increase  in  GNR  width.  This  increase  can  be 
attributed  to  the  larger  area  of  overlap  between  the  crossed-over  GNRs.  More  importantly,  the 
figure  also  shows  that  for  all  cases,  the  total  thermal  energy  exchange  (and  thus  /\ ,)  increases  as 
well  with  GNR  length  for  a  constant  GNR  width.  As  discussed  later,  this  increase  is  attributed 
to  an  increased  role  of  GNR  edges  towards  inter-graphene  coupling  with  the  increase  in  length. 
It  is  also  evident  from  this  figure  that  A\  has  a  sub-Vi  power  law  like  scaling  with  respect  to  the 
GNR  length. 

Figure  2(b)  plots  the  discussed  data  from  the  perspective  of  an  area-normalized  thermal 
conductance  (A2).  Interestingly,  this  figure  details  several  interesting  and  key  observations.  For 
physically  interacting  GNRs  as  investigated  in  this  study,  one  would  expect  the  normalized 
conductance  to  be  somewhat  independent  of  the  overlapping  area  and  nano-ribbons  parameters. 
However,  the  figure  clearly  shows  that  a  single  value  of  normalized  thermal  conductance  is  not 
valid  but  ranges  between  ~20-70  MW/m2  K  for  the  studied  cases.  Moreover,  the  figure  shows 
that  the  A2  increases  with  GNR  length  (for  a  constant  width)  and  decreases  with  GNR  width  (for 
a  constant  length).  The  decrease  in  A2  for  a  constant  GNR  length  suggests  that  if  one  increases 
the  overlap  area  by  a  constant  factor  (by  increasing  GNR  width),  the  increase  in  thermal  energy 
exchange  is  less  than  proportional  to  the  contact  area  increase  (Please  see  below  for  further 
discussion).  Overall,  it  is  quite  clear  from  figure  2  that  cross-plane  thermal  conductance  values 
(both  Ai  and  A2)  between  GNRs  are  strongly  dependent  on  their  dimensions. 

Thermal  conductance  versus  aspect  ratio 

From  the  observed  Ai  and  A2  dependences  on  L  and  W  of  the  GNRs  in  contact,  it  is  reasonable 
to  propose  that  the  aspect  ratio  {rf)  should  be  an  important  parameter  associated  with  thermal 
transport  at  the  contact  interface.  In  this  context,  thermal  conductance  data  in  figure  2  was  re¬ 
analyzed  as  a  function  of  aspect  ratio  and  is  shown  in  figure  3  for  area-normalized  thermal 
conductance.  A  similar  plot  for  total  thermal  conductance  versus  //  did  not  add  any  additional 
insights  and  is  provided  in  the  supporting  information  S2  for  reader’s  interest.  The  aspect  ratio 
for  the  studied  cases  was  calculated  using  the  following  equation 
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n  = 


L' 


A 


GNR 


(3) 


where  L'  is  the  length  (in  nm)  of  the  GNR  excluding  fixed  ends  and  AGNR  is  the  effective  area 
(in  nm2)  of  contact  between  GNRs  (Please  see  supporting  information  SI  for  calculation 
°f  ^gnr)- 

It  is  fascinating  to  observe  that  when  A2  is  plotted  with  respect  to  aspect  ratio  (unlike 
length  or  width),  all  the  data  points  fall  on  a  single  master  plot  up  to  the  aspect  ratio  of  ~20.  The 
figure  clearly  shows  that  the  aspect  ratio  of  the  GNR  is,  in  fact,  the  governing  parameter  which 
determines  the  extent  of  thermal  energy  exchange  between  GNRs.  For  example,  in  accordance 
with  the  trend  in  figure  3,  for  two  different  size  GNRs  with  same  aspect  ratio  (e.g.  with 
dimensions  WxL  =  4x40n m  and  WxL  =  8  x  80 nm),  we  found  A2  to  be  45.5  and  44.5  MW/ 
m2K  respectively.  In  addition,  the  associated  per  unit  area  thermal  energy  exchange  was 
observed  to  be  0.50  and  0.53  MW/m2  respectively. 

Figure  3  also  shows  a  power  law  scaling  guideline  of  Vi.  While  at  low  aspect  ratios,  the 
scaling  of  Vi  seems  fitting,  it  shows  a  considerable  deviation  at  higher  aspect  ratios.  In  order  to 
get  a  better  understanding  of  normalized  conductance  with  respect  to  aspect  ratio,  we  tested  the 
fit  of  the  data  with  several  functional  forms.  The  two  most  appropriate  fits  among  those  are 
shown  in  figures  4(a)  and  4(b),  namely,  the  power  fitting  and  a  three-parameter  logarithmic 
fitting,  respectively.  The  corresponding  equations  and  fitting  exponents  are  also  shown  in  the 
respective  figures.  For  power  fitting  (figure  4(a)),  a  much  smaller  exponent  (0.38)  than  Vi  is 
observed  when  fitted  against  all  data  points.  Moreover,  it  is  easy  to  observe  clear  deviations  of 
the  data  points  with  respect  to  the  power-law  fitting  line.  Most  of  the  middle  data  points  reside 
above  the  fitting  line,  while  the  ones  on  the  extremes  are  below  the  line.  On  the  other  hand,  the 
logarithmic  fitting  (figure  4(b))  seems  to  provide  a  better  overall  fitting.  The  quality  of  the  fit 
can  also  be  gauged  by  comparing  the  fitting  exponents  for  a  full  dataset  fit  with  fittings  done 
over  a  smaller  window.  Such  fittings  are  reported  in  supporting  information  S3  for  power  and 
logarithmic  fittings,  respectively.  It  is  seen  that  the  logarithmic  fitting  on  partial  data  fits  results 
in  the  same  coefficients  when  compared  to  a  full  fit  while  a  clear  distinction  is  observed  for 
power  law  fitting  (coefficient  a  monotonically  increases  from  16.13  to  21.89  while  coefficient  b 
decreases  from  0.5  to  0.31),  further  validating  the  previous  discussion.  Moreover,  based  on 
functional  form,  it  is  easy  to  comprehend  that  a  logarithmic  fitting  leads  to  a  better  asymptotic 
fit  at  larger  aspect  ratios  (~20-100),  where  the  thermal  conductance  is  expected  to  saturate. 


Proposed  phenomenological  model 

We  propose  the  following  phenomenological  model  to  qualitatively  appreciate  the  dependence 
of  ITC  on  aspect  ratio  as  well  as  its  logarithmic  nature.  In  our  model,  the  ITC  ( A2 )  is  effectively 
a  weighted  average  of  two  limiting  thermal  conductance  values,  4,  and  A0.  Here,  sub-scripts  i 
and  o  refer  to  the  thermal  conductance  due  to  inner  and  outer  interacting  regions,  respectively, 
and  are  also  shown  in  figure  5.  In  our  model,  A,  refers  to  the  intrinsic  thermal  conductance 
across  bi-layer  graphene  in  the  limit  of  no  edges,  i.e.  rj  ~  1  or  large  bi-layer  graphene  flakes.  On 
the  other  hand,  A0  corresponds  to  thermal  conductance,  when  a  GNR  ‘edge’  dominates  the 
thermal  interaction  with  the  other  GNR.  Due  to  the  greater  flexibility,  i.e.  less  constrained  out- 
of-plane  vibrations  of  edge  atoms  than  atoms  in  the  inner  interaction  region,  we  assume  that 
edges  can  couple  with  other  GNR  to  a  stronger  degree,  resulting  in  a  higher  limiting  value  of  A0 
than  Al  The  differences  in  vibrational  properties  of  inner  and  edge  atoms  are  also  highlighted  in 
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Figure  4.  Graphs  showing  (a)  power  fitting  and  (b)  logarithmic  fitting  of  area- 
normalized  thermal  conductance  (A2)  data  sets  with  respect  to  aspect  ratio  0/).  The 
equations  along  with  derived  parameters  are  also  listed  in  the  respective  graphs.  Please 
refer  to  figures  S3a  and  S3b  in  the  supporting  information  for  comparing  the  derived 
parameters  with  partial  fitting  of  the  data  set. 


figure  6  as  a  power  spectra  for  a  representative  system,  calculated  from  the  Fourier  transform  of 
the  velocity  auto-correlation  function.  The  figure  clearly  shows  the  shift  of  vibrational  spectrum 
to  lower  frequencies  for  edge  atoms  (especially  ZO  mode  which  occurs  around  ~20  THz  and  is 
responsible  for  out  of  plane  vibrations),  validating  the  greater  degree  of  vibrational  flexibility. 

We  also  foresee  a  length  ld  (shown  in  figure  5)  from  the  GNR  edge  over  which  the  thermal 
energy  exchange  is  dominated  through  edges,  i.e.  A0.  We  will  refer  to  ld  as  decay  length  as  it 
indicates  a  length  scale  over  which  the  effect  of  edge  vibrational  characteristics  is  noticeable. 
Within  this  framework,  the  effective  ITC  (A2)  can  be  calculated  using  the  following  equation  of 
geometric  averaging  in  terms  of  ld,  W,  Aj  and  A0. 
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Inner  Region  (A) 


Outer  Region  (A0) 


Figure  5.  Schematic  showing  interaction  (or  overlap)  region  and  identifying  the  inner 
and  outer  regions  and  their  respective  thermal  conductance  parameters  (Ax  and  A0).  Zd  is 
the  decay  length  over  which  the  effect  of  the  edges  become  negligible. 


Frequency  (THz) 


Figure  6.  Plot  of  phonon  density  of  states  versus  frequency  for  carbon  atoms  in  inner 
(black)  and  outer  (red)  interaction  regions  for  a  representative  system  (L  =  40  nm  and 
W=  10  nm)  with  aspect  ratio  of  4. 


^2 


A\  x  ( W-2ldf  +  2 A0(Wxld)  +  2A0((W-2ld)  x  /d) 

W2 


In  the  limit,  when  /d  is  significantly  smaller  than  W,  equation  (4)  can  be  simplified  as 
A2  =  Ai  +  ^(2 A0  — A;). 


(4) 

(5) 


We  assume  that  the  decay  length  is  not  constant  and  is  a  function  of  both  L  and  W.  It  is 
expected  that  for  a  constant  W,  ld  oc/(L)  as  the  effect  of  edges  towards  thermal  energy  exchange 
is  expected  to  penetrate  deeper  into  the  interacting  region  with  increasing  distance  (L)  between 
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the  interaction  region  and  the  fixed  (constrained)  GNR  ends,  due  to  the  greater  degree  of  edges’ 
out-of-plane  flexibility.  Moreover,  when  //  becomes  unity,  ld  is  expected  to  vanish  since  there 
would  not  be  edges,  and  hence  no  decay  length.  One  function  which  satisfies  both  conditions  is 
a  logarithmic  function.  Here,  ld  oc  In (L/W)  or  ld  oc  ln(//).  This  functional  form  suggests  that  for 
constant  W,  the  decay  length  increases  ‘logarithmically’  with  distance  between  the  fixed 
constraints. 

From  the  ITC  results,  it  is  also  observed  that  the  relative  contributions  from  A;  and  A0 
remains  similar  at  certain  constant  rj.  In  order  to  accommodate  that,  ld  should  linearly  increase 
with  W  when  aspect  ratio  is  constant.  This  linear  proportionality  of  ld  with  W  is  also  in 
accordance  with  the  need  of  additional  length  dimension  in  ld  oc  with  respect  to 

dimensional  equality  of  left-  and  right-hand  side  of  the  equations.  Overall,  we  propose  that  ld  oc 
Wln(ti )  or  ldIW=iAn{rf)  where  yu  is  a  dimensionless  coefficient  of  proportionality.  Putting  this 
form  in  equation  (5),  we  get 

A2  =  +  2/u(2A0  — A/)  In  (rf).  (6) 

As  can  be  seen,  equation  (6)  is  very  similar  to  fitted  equation  in  figure  4(b)  and  showcases 
the  observed  logarithmic  dependence  when  ld  «  W  and  secondary  terms  can  be  ignored,  which 
in  general  holds  true  for  small  aspect  ratio  systems.  When,  ld  is  comparable  to  W  (higher  aspect 
ratios),  the  secondary  terms  originating  from  4  comers  (see  equation  (4))  should  not  be  ignored. 
In  such  cases,  equation  (6)  can  be  derived  to  be 

A2  =  Ai  +  2fi ( 2A0  -Ai)Mri)  -4 M2(A0  -A^ln^))2  .  (7) 

Although  the  functional  form  in  equation  (7)  is  not  similar  to  the  fitted  form  in  figure  4(b), 
it  does  provide  the  saturation  of  A2  at  higher  aspect  ratios.  With  regards  to  prediction  of  A;  in 
limiting  case  of  tj  ~  1  (negligible  contributions  from  the  edges),  the  extrapolation  of  three- 
parameter  logarithmic  fitting  lead  to  estimated  value  of  15.71  MW/m2-K,  which  is  also  in  good 
agreement  with  what  is  reported  in  the  literature  [28].  Furthermore,  based  on  observed  trends, 
we  expect  the  intrinsic  value  of  ITC  A0  to  be  70  MW/m2-K. 

We  should  point  out  that  the  derivation  of  our  proposed  model  is  phenomenological  in 
nature  and  is  proposed  for  appreciating  observed  trends  and  provide  a  qualitative  understanding 
of  the  features  which  are  important  across  the  interacting  zone. 

Effect  of  edge  roughness 

Previous  discussion  suggests  that  the  vibrational  freedom  of  GNR  edges  do  play  an  important 
role  in  determining  the  decay  length  and  eventually  the  total  thermal  energy  exchange  across 
the  interacting  bi-layer  GNR  interface.  Hence,  it  is  imperative  to  investigate  the  potential  role  of 
edge  roughness  towards  out-of-plane  thermal  energy  exchange.  In  this  regard,  we  simulated 
ITC  across  a  representative  rough  GNR  for  different  lengths  L  between  fixed  constraints,  but 
similar  interacting  area,  as  depicted  in  figure  7.  As  seen  in  the  schematic,  although  the  GNR  is 
significantly  rough  with  regards  to  transmitting  in-plane  phonons,  the  interacting  zone  was  kept 
to  be  pristine  and  crossing  at  90°  in  order  to  compare  the  results  with  rectangular  (non-rough) 
GNR  studies.  Figure  7  also  shows  the  calculated  thermal  conductance  (A2)  values  across  rough 
GNRs  as  a  function  of  aspect  ratio  and  compare  them  with  A2  values  corresponding  to 
rectangular  GNRs  (shown  as  partially  transparent  data  points).  The  figure  shows  that  the  rough 
GNR  thermal  conductance  values  match  excellently  with  those  with  rectangular  GNRs, 
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Figure  7.  Plot  of  area-normalized  interface  thermal  conductance  as  a  function  of  aspect 
ratio  (pink  hexagons),  showing  the  effect  of  overall  graphene  roughness  on  the  thermal 
energy  exchange  across  the  interface.  The  top-inset  shows  both  the  overall  structure  of 
rough  GNR  as  well  as  the  structure  interaction  region.  The  partially  transparent  data 
points  correspond  to  thermal  conductance  values  of  studied  rectangular  GNRs  as 
discussed  in  figure  3. 

suggesting  that  overall  zigzag  morphology  of  the  GNR  do  not  play  a  role  towards  out  of  plane 
thermal  conductance  across  GNRs.  In  other  words,  a  decrease  in  overall  in-plane  thermal 
conduction  properties  (for  example,  due  to  roughness)  does  not  directly  imply  better  out-of- 
plane  conduction.  Furthermore,  it  is  observed  that  the  interacting  zone  between  GNRs  is  the 
same  for  both  rough  and  rectangular  GNRs  and  leads  to  similar  values  of  A2.  Based  on  these 
observations,  we  argue  that  the  out-of-plane  thermal  energy  exchange  depends  predominantly 
on  local  vibrational  characteristics  of  the  interacting  atoms,  rather  than  overall  roughness 
features  of  GNRs.  Comparing  rough  and  rectangular  GNR  cases,  it  can  also  be  argued  that  the 
said  local  vibrational  characteristics  depends  more  significantly  on  the  L  (distance  between  the 
constraints)  and  W  (GNR  width),  but  not  so  much  on  overall  roughness  features  of  GNRs. 
However,  we  hypothesize  that  if  the  interaction  region  is  modeled  to  be  rough  as  well  (which  is 
not  investigated  here),  a  different  scaling  could  emerge  as  it  would  modify  the  decay  length 
characteristics,  such  as  switching  from  a  pristine  one-dimensional  character  (perpendicular  to 
the  edge)  to  some  different  decay  characteristics  with  higher  degree  of  complexity  (such  as  its 
dependence  on  degree  of  roughness). 


Effect  of  temperature 

To  investigate  the  effect  of  temperature  on  thermal  conductance  across  the  GNR  contact,  we 
simulated  a  representative  system  with  L  and  W  correspond  to  60  nm  and  lOnm,  respectively  at 
different  temperatures  values  for  hot  and  cold  thermostats,  using  the  same  NEMD  methodology 
discussed  in  the  text.  In  this  context,  figure  8  plots  the  calculated  ITC  (A2)  values  for  different 
equilibrated  temperature.  It  is  clear  from  the  figure  that  the  A2  increases  with  temperature  and 
could  be  explained  using  an  argument  based  on  the  number  of  interacting  phonons.  While  the 
aspect  ratio  was  not  altered  for  this  study  (effect  of  temperature),  the  number  of  phonons 
interacting  across  the  interface  increase  (decrease)  with  increase  (decrease)  in  temperature. 
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Figure  8.  Plot  of  area-normalized  interface  thermal  conductance  versus  temperature  as 
studied  for  a  representative  cross  shaped  geometry  of  60  nm  long  and  10  nm 
wide  GNRs. 


z  Normalized  Strain  Ratio 

Figure  9.  Plot  of  area-normalized  thermal  conductance  across  40  nm  long  and  10  nm 
wide  GNR  as  a  function  of  compressive  (<1.0)  and  tensile  (>1.0)  strain  values. 

Furthermore,  the  mean  free  path  of  the  phonons  becomes  longer  at  a  lower  temperature  due  to 
lower  number  of  phonons  and  a  lesser  number  of  in-plane  anharmonic  phonon-phonon 
scattering  events.  We  foresee  that  such  an  increase  in  the  mean  free  path  would  lower  the 
phonon  interaction  time  across  the  GNR  contact  and  vice-versa.  We  believe  that  both  discussed 
points  contribute  towards  increasing  (decreasing)  the  intrinsic  values  of  both  and  A0  with 
increase  (decrease)  in  temperature,  which  in  turn  results  in  shown  variation  (figure  8)  in 
effective  ITC  with  temperature. 


Effect  of  tensile  &  compressive  strain 

Figure  9  shows  the  predicted  values  of  normalized  thermal  conductance  (A2)  for  different 
degrees  of  compressive  and  tensile  strain  with  respect  to  unstrained  GNR.  Also  shown  in  the 
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Figure  10.  Plot  of  phonon  density  of  states  versus  frequency  for  carbon  atoms  for  a 
representative  system  (L  =  40nm  and  W=  10  nm)  with  aspect  ratio  of  4,  focusing  on 
shifting  of  ZO  vibrational  mode  peak  for  different  degree  of  applied  tensile  strain. 

figure  are  a  few  schematic  snapshots  of  the  strained  systems  in  the  equilibrated  geometry  for  a 
clearer  visualization.  The  figure  shows  two  trends  which  are  discussed  below. 

Under  tensile  strain,  a  significant  decrease  in  A2  was  observed.  Here,  several  factors  play  a 
role  which  could  potentially  alter  the  observed  values.  First,  the  areal  density  of  interacting 
carbon  atoms  decreases  with  tensile  strain.  For  example,  with  a  relative  strain  of  1.1,  the  areal 
carbon  atom  density  decreases  by  a  factor  of  10%  on  both  GNRs,  which  in  turn  results  in  a 
~20%  decrease  in  interaction  energy  for  the  same  overlap  area.  Secondly,  as  discussed  above, 
the  increase  in  aspect  ratio  (~1.1  times)  is  also  expected  to  change  (~5%  increase  from  figure  3 
data)  the  thermal  conductance.  While  these  combined  factors  could  lead  to  a  predicted 
~  10-20%  decrease  in  thermal  conductance  values,  the  observed  five-fold  decrease  can  be 
explained  using  a  bond  stiffening  argument.  It  is  expected  that  GNRs  become  notably  stiffer 
with  tensile  strain  due  to  stiffening  of  C-C  bonds.  This  effect  of  bond  stiffening  is  clearly  shown 
in  figure  10  for  different  applied  tensile  strain  values  for  the  representative  studied  system 
(L-  40 nm  and  W=  \  0  nm),  focusing  on  out-of-plane  optical  mode  (ZO)  [50],  occurring 
~20THz.  With  increasing  values  of  in-plane  strain,  the  figure  clearly  shows  a  shift  of 
characteristic  vibrational  peak  towards  higher  frequencies.  This  apparent  shift  to  the  higher 
frequency  (higher  energy)  leads  to  a  lower  number  of  interacting  phonons  (lower  amplitude)  at 
studied  temperature.  Furthermore,  the  stiffening  is  expected  to  increase  the  group  velocity  of  the 
phonons  and  increase  their  effective  mean  free  path.  This,  in  turn,  relatively  lessens  the  number 
of  anharmonic  phonon-phonon  scattering  events,  also  lowering  the  effective  phonon  interaction 
time  across  the  junction.  We  argue  that  both  factors  (lower  number  of  phonons  and  lesser 
interaction  time  due  to  bond  stiffening)  contribute  towards  observed  lowering  of  thermal 
conductance  with  increasing  tensile  strain. 

For  applied  compressive  strain,  it  was  observed  that  the  conductance  was  not  noticeably 
affected.  When  compressive  strain  was  used  to  shorten  the  GNRs  along  its  length  direction  (by 
shortening  all  bond  lengths  at  beginning  of  the  strain),  we  observed  a  subsequent  relaxation  of 
the  GNRs  by  a  resetting  of  the  bonds  to  their  original  unstrained  lengths.  To  accommodate  such 
relaxation,  the  nano-ribbons  buckled  to  some  degree  (see  schematic  in  figure  9  corresponding  to 
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a  normalized  strain  ratio  of  0.94).  Additionally,  the  contact  region  was  observed  to  be  relatively 
flat  for  all  buckled  cases.  Since  the  areal  density  of  the  carbon  atoms  and  the  contact  area 
remained  very  similar  even  after  bucking,  we  observed  that  the  thermal  energy  exchange  was 
not  noticeably  affected,  resulting  in  very  similar  values  of  ITC  (A2).  The  observed  invariance  in 
A2  is  indicative  of  similar  contributions  from  and  A0  for  different  degrees  of  partial  GNR 
buckling.  Similar  to  indifferences  in  A2  due  to  studied  roughness,  the  simulations  further 
substantiate  that  out-of-plane  thermal  conductance  predominantly  depends  upon  phonon 
characteristics  of  the  interacting  region  and  is  indifferent  to  overall  GNR  conformations  such  as 
partial  GNR  buckling  and  their  roughness  characteristics  as  discussed  before. 


Summary  and  conclusions 

The  observed  thermal  transport  dependences  in  GNR  crossed-over  contacts  from  MD 
simulations  demonstrate  a  few  important  points.  GNR  aspect  ratio  is  found  to  be  a  key  physical 
characteristic  which  governs  the  magnitude  of  thermal  energy  exchange  and  thus,  thermal 
conductance.  The  dependence  of  aspect  ratio,  found  to  be  logarithmic  in  nature,  is  derived  using 
a  phenomenological  model,  decoupling  the  contributions  from  GNR  edges  and  inner  atoms 
within  the  interacting  zone.  The  invariance  of  thermal  conductance  with  respect  to  studied 
rough  GNR  indicates  that  the  out-of-plane  conductance  (A2)  is  determined  by  the  edge 
characteristics  within  the  interacting  zone  and  does  not  noticeably  depend  on  overall  GNR 
roughness  features.  Furthermore,  the  temperature  and  tensile  strain  dependence  on  calculated 
values  of  thermal  conductance  is  shown  to  be  governed  by  the  number  of  interacting  phonons  as 
well  as  their  retention  time  across  the  interacting  zone. 

We  should  duly  note  that  because  of  the  classical  nature  of  the  performed  MD  simulations, 
the  effect  of  electron-phonon  interactions  was  not  included  in  this  study.  It  is  likely  that  for 
GNRs  with  metallic  or  small  band  gap  semiconducting  character,  the  interfacial  thermal 
conductance  will  be  strongly  influenced  by  the  electronic  conduction  as  well  as 
electron-phonon  interactions  in  the  contact.  Furthermore,  it  is  important  to  note  that  as  the 
conductance  values  were  calculated  in  a  suspended  GNR-pair  geometry,  their  direct 
incorporation  into  design  rules  for  supported-GNR  based  nano-electronics  should  be  performed 
with  cautious,  as  graphene-substrate  interactions  (which  are  neglected  here)  should  also  play  an 
important  role  towards  determining  the  extent  of  thermal  energy  exchange  and  temperature 
distribution  within  the  system.  Nevertheless,  we  believe  that  our  study  provides  the  first  few  key 
steps  in  that  direction  from  the  perspective  of  which  intrinsic  GNR  parameters  are  of  primary 
importance  and  the  extent  of  modulation  in  ITC  across  them. 
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